A Multiphase Coupled Hydrodynamic Model for Fate and Transport Simulation of Polycyclic Aromatic Hydrocarbons in a Semi-Closed Narrow Bay

With their unique geographical characteristics, semi-closed narrow bays are important places for human survival but vulnerable to pollution. Because pollutants (polycyclic aromatic hydrocarbons, PAHs) migrate and undergo transformation through a dynamic mechanism in bays of this type, environmental authorities have formulated a series of effective measures for pollution prevention and control, but these are difficult to realize. Based on monitoring and historical data, a multiphase-coupled hydrodynamic model combined with a carcinogenic risk-assessment model was able to solve the challenging environmental problem. Results showed that the hydrodynamic condition in the semi-closed narrow bay was very complex. A weaker hydrodynamic force had an adverse influence on the diffusion of pollutants, further amplified in part by the head of the semi-closed narrow bay, resulting in a higher ecological risk. The prediction results indicated that the total amount of PAHs transported from seawater to sediments was about 4.7 × 1013 ng/year, which might cause serious threats to aquaculture or human health.


Introduction
With their unique geographical environment, superior natural conditions, and rich marine resources, bays are the best sites for human development and utilization, especially the semiclosed narrow bays [1]. Because of their limited width and extended length, they are ideal places for aquaculture, transport and harbor industries [2]. They play an important role in human survival and developing and utilizing marine resources [3]. However, semi-closed narrow bays also have negative properties [4]. Located in the sea-land interaction area, the substance, energy and physicochemical characteristics of the two systems are integrated closely [5], which makes narrow semi-closed bays very vulnerable to the influence of massive pollutants from the land [6]. The features of narrowness and length cause difficulties in self-purification and diffusion of pollutants [7], creating weak hydrodynamic conditions. Once pollutants are discharged, they will cause long-term pollution which will pose a serious threat to the health and stability of the marine ecosystem [8].
In recent years, the use of fossil fuels and the landfill and incineration of domestic waste have constantly introduced persistent toxic substances or PAHs into the environment in China's coastal areas [9][10][11]. With the strong teratogenic, carcinogenic, and lethal effects of toxins on the larvae and embryos of algae, shrimp and shellfish [12][13][14], PAHs could seriously endanger the growth, development and reproduction of aquatic organisms [15]. In addition, PAHs are not easy to metabolize and decompose, resulting in biological enrichment through the food chain, which would cause harm to fish, marine mammals and even human beings [3,16]. Consequently, the pollution situation in semi-closed, narrow

Construction of the Model
With the help of Delft3D, a multiphase coupled hydrodynamic model of the PAHs in Pulandian Bay was constructed. The dynamic equation was solved by using the explicit and implicit alternating numerical integration method to simulate the migration and transformation process of PAHs in the marine environment [29][30][31]. Figure 2 is the schematic diagram of the model.
∂V ∂t PAHs Kinetic diffusion: Exchange dynamics of PAHs between seawater and atmosphere: k vol = 1/( 1 k l + R g (Tem + 273.15) N g Pk g e (a 1 +a 2 /(Tem+273.15)) ) (7) Exchange dynamics of PAHs between seawater and seabed: Comprehensive decomposition (degradation) process: Stokes formula (deposition): where a 1 was the temperature coefficient for volatization entropy, a 2 was the temperature coefficient for volatilization enthalpy, C D was the coefficient of wind-stress, C i , C equ , C sed , C w , C par , and C d were the total concentration, equilibrium concentration, sediment concentration, concentration in upper overlying water, particle adsorption state concentration, and free dissolved state concentration of PAHs, d was the size of suspended particles, D sw was the diffusion coefficient, f was the Coriolis parameter, g was the acceleration of gravity, H was the water depth, K 0 was the zero-order reaction constant, K 1 was the first-order reaction constant, K t was the thermodynamic constant, k l was the transfer coefficient for the liquid film, k g was the transfer coefficient for the gas film, k vol was the seawater-air exchange coefficient, L sw was the thickness of the overlying water layer, Ng was the number of moles in a m 3 gas, P was the atmospheric pressure, Q was the pollutant discharge intensity, R was the decomposition rate of PAHs, R D was the diffusion flux, R g was the gas constant, R vol was the seawater-air exchange rate, R p was the permeation flux, t was time, Tem was the ambient temperature, U,V were the velocity components in the X,Y directions, V 0 was the PAHs sedimentation rate, V e was the seawater viscosity coefficient, v p seepage velocity, W was the wind-speed at 10m above water surface, σtx and σty were the turbulent diffusion coefficient in X and Y direction, ζ was the water level, ψ was the Chézy coefficient, τ (including τ x , τ y ) was wind-stress, ϕ was the sediment porosity, ω was the concentration of adsorbed PAHs, ρ 0 was the density of sea water, and ρ T was the density of suspended particles.

Toxicity Analysis Model of PAHs
By ingestion and dermal uptake, PAHs are able to enter the body and spread to various organs, such as the brain, muscles, fat, liver, and kidneys, through alveoli, the digestive system, and blood [17,26]. Therefore, the US Environmental Protection Agency method [32] was used to estimate incremental lifetime cancer risks (ILCRs) to assess potential health risks to humans posed by PAHs.
The ILCRs for the ingestion and dermal uptake were calculated using the equations above [33,34].
Considering that PAHs are a mixture with many components, many physical or chemical differences exist. To ensure the consistency of PAHs simulation results, the toxicityequivalent-concentration method was used in this study [35]. For toxicity-equivalent factors of 16 PAHs, see Table 1. The total concentration of 16 PAHs was replaced by a Benzo(a)pyrene toxicity-equivalent concentration.
where C i was the concentration of individual species of PAHs in a calculation grid (in Figure 3), TEF i was the toxicity-equivalent factor of individual species of PAHs (Table 1.), ΣPAHs was the total toxicity equivalent of 16 PAHs relative to BAP (ng/g). ABS was the dermal adsorption fraction, AF was the dermal adherence factor (mg/(cm 2 ·h)), AT was the average life span (y), BW was the average body weight (kg), C S was the ΣPAHs concentration (ng/g) in sediment (i.e., ΣPAHs), C SF was the carcinogenic slope factor for ingestion or dermal uptake (mg/(kg·d)), ED was exposure duration (y), EF was exposure frequency (d/y), IR Ingestion was the intake rate (mg/d), and SA was the dermal exposure area (cm 2 ).  The ILCR method mentioned above was used to evaluate the health risk caused by PAHs. The ILCR parameters used in the equation are shown in Tables 1 and 2. ∑ ILCR = ILCR ingestion + ILCR dermal (16)  The ILCR value was divided into three levels, according to the U.S. Environmental Protection Agency's guidelines on health risk assessment of PAHs. When the ILCR value < 10 -6 , it indicates that the health risk is low or negligible; when the ILCR value is between 10 -6 -10 -4 , it indicates a potential carcinogenic risk; when the ILCR value >10 -4 , it indicates a high carcinogenic risk.

Model Parameter Setting
The boundary of the study area was extracted from the satellite remote sensing image. ArcGIS (10.3, ESRI) was applied to vectorize the terrain, water depth, and boundary and determine the open and closed boundary as well as the calculation domain. The dynamic process of pollutants in a narrow and semi-closed bay was so complex that orthogonal curved grids were used, which could be closer to the terrain boundary. To cover every small corner of the complex topography of Pulandian Bay, the whole calculation domain was compartmentalized into 32,332 (274 × 118) orthogonal curved grids ( Figure 3).
The water depth in Pulandian Bay was shallow, ranging from 5 to 25 m, and the topographic relief was smooth and level. A two-dimensional ocean hydrodynamics model was suitable. The model-time step was set as 5.0 min. Gravity, geostrophic deflection force, and tidal effect were the main driving forces of the model. The gravity acceleration g was 9.81 m/s 2 , the seawater density was 1024.2 kg/m 3 , the diffusion and viscosity coefficients were 10.1m 2 /s and 1.0 × 10 −6 m 2 /s. Tidal harmonic constants were obtained from the harmonic analysis based on 2015 tidal data (including, K 1 , M 2 , O 1 , S 2 , P 1 , Q 1 , M 1 , N 2 , K 2 ). The primary parameters of the model are shown in Table 3.

Initial Conditions
The dynamic process of pollutants in the narrow and semi-closed bay was very complex. The hydrodynamics force, the physical and chemical reactions of air-seawater, seawater-sediment, and the degradation of pollutants themselves lead to large fluctuations in their equilibrium state in different media (such as air, water, or sediment). To ensure the accuracy of the numerical calculation model simulation, the calculation scheme of "hot start" was implemented instead of "cold start." Therefore, the initial concentration of PAHs in Pulandian Bay was set as

Initial Conditions
The dynamic process of pollutants in the narrow and semi-closed bay was very complex. The hydrodynamics force, the physical and chemical reactions of air-seawater, seawater-sediment, and the degradation of pollutants themselves lead to large fluctuations in their equilibrium state in different media (such as air, water, or sediment). To ensure the accuracy of the numerical calculation model simulation, the calculation scheme of "hot start" was implemented instead of "cold start." Therefore, the initial concentration of PAHs in Pulandian Bay was set as ϛi and ξ (according to the monitoring results of PAHs).
The pollution source intensity in the model showed the equivalent average value of 120 μg/L·(s −1 m 3 ) [28]. The position of the discharge sources is shown in Figure 1. Once the simulation started, a series of migration and transformation processes of PAHs, such as atmosphere-seawater exchange, seawater-sediment exchange, and sedimentation and degradation began. Synchronously, the mixing process of seawater both inside and outside the bay occurred, which was the dominant function for contaminant diffusion.

Accuracy Verification of the Model
To ensure the multiphase coupled hydrodynamic model could effectively simulate the fate and transport of PAHs in Pulandian Bay, accuracy verification was critical. Due to the massive amount and complexity of calculation, not only a strong mathematical and computer programming ability was required, but also a large number of long-term historical data needed to be used for calibrating the model parameters and improving its accuracy. In this study, three aspects were taken into consideration for the accuracy verification of the model. The first was to ensure the accuracy of the on-site monitoring data of PAHs and hydrology; the second was to calibrate the model parameters during the process of model construction; the third was to verify the accuracy of the simulation results after the calculation.

Collection and Experimental Analysis
In the summer of 2020, the investigation and sampling of hydrology, water quality and sediment in Pulandian Bay through a field survey were carried out. A total of 32 monitoring stations were set up. To verify the accuracy of the model simulation, three water-quality monitoring stations, three hydrological-monitoring stations (S3, S18 and S20), and two tide-level observation stations (S11 and S15) were set up. The layout of the sampling stations is shown in Figure 1. A grab sampler (iron) was used to get the surface sediment (0-10 cm) samples, and a layered water sampler (stainless steel, 5 L) was used to get seawater samples. Hydrologic observations (flow velocity and direction) were collected every 10 min by using an Acoustic Doppler Current Profiler (Model: WHS600, depth unit size: 600 kHz).
Seawater samples were enriched by an envi-C18 solid phase membrane disk (flow rate: 12 mL/min−30 mL/min, vacuum filtration) extraction, wrapped in aluminum foil, sealed, and stored before testing (SPE samples). Sediment samples were freeze-dried (vacuum freeze dryer, 24 h) and ground (180 μm) after removal of impurities (e.g., plant residues and animal). Brown glass bottles were used for storage and then brought back for analysis.

Initial Conditions
The dynamic process of pollutants in the narrow and semi-closed bay was very complex. The hydrodynamics force, the physical and chemical reactions of air-seawater, seawater-sediment, and the degradation of pollutants themselves lead to large fluctuations in their equilibrium state in different media (such as air, water, or sediment). To ensure the accuracy of the numerical calculation model simulation, the calculation scheme of "hot start" was implemented instead of "cold start." Therefore, the initial concentration of PAHs in Pulandian Bay was set as ϛi and ξ (according to the monitoring results of PAHs).
The pollution source intensity in the model showed the equivalent average value of 120 μg/L·(s −1 m 3 ) [28]. The position of the discharge sources is shown in Figure 1. Once the simulation started, a series of migration and transformation processes of PAHs, such as atmosphere-seawater exchange, seawater-sediment exchange, and sedimentation and degradation began. Synchronously, the mixing process of seawater both inside and outside the bay occurred, which was the dominant function for contaminant diffusion.

Accuracy Verification of the Model
To ensure the multiphase coupled hydrodynamic model could effectively simulate the fate and transport of PAHs in Pulandian Bay, accuracy verification was critical. Due to the massive amount and complexity of calculation, not only a strong mathematical and computer programming ability was required, but also a large number of long-term historical data needed to be used for calibrating the model parameters and improving its accuracy. In this study, three aspects were taken into consideration for the accuracy verification of the model. The first was to ensure the accuracy of the on-site monitoring data of PAHs and hydrology; the second was to calibrate the model parameters during the process of model construction; the third was to verify the accuracy of the simulation results after the calculation.

Collection and Experimental Analysis
In the summer of 2020, the investigation and sampling of hydrology, water quality and sediment in Pulandian Bay through a field survey were carried out. A total of 32 monitoring stations were set up. To verify the accuracy of the model simulation, three water-quality monitoring stations, three hydrological-monitoring stations (S3, S18 and S20), and two tide-level observation stations (S11 and S15) were set up. The layout of the sampling stations is shown in Figure 1. A grab sampler (iron) was used to get the surface sediment (0-10 cm) samples, and a layered water sampler (stainless steel, 5 L) was used to get seawater samples. Hydrologic observations (flow velocity and direction) were collected every 10 min by using an Acoustic Doppler Current Profiler (Model: WHS600, depth unit size: 600 kHz).
Seawater samples were enriched by an envi-C18 solid phase membrane disk (flow rate: 12 mL/min−30 mL/min, vacuum filtration) extraction, wrapped in aluminum foil, sealed, and stored before testing (SPE samples). Sediment samples were freeze-dried (vacuum freeze dryer, 24 h) and ground (180 μm) after removal of impurities (e.g., plant residues and animal). Brown glass bottles were used for storage and then brought back for analysis.
The pollution source intensity in the model showed the equivalent average value of 120 µg/L·(s −1 m 3 ) [28]. The position of the discharge sources is shown in Figure 1. Once the simulation started, a series of migration and transformation processes of PAHs, such as atmosphere-seawater exchange, seawater-sediment exchange, and sedimentation and degradation began. Synchronously, the mixing process of seawater both inside and outside the bay occurred, which was the dominant function for contaminant diffusion.

Accuracy Verification of the Model
To ensure the multiphase coupled hydrodynamic model could effectively simulate the fate and transport of PAHs in Pulandian Bay, accuracy verification was critical. Due to the massive amount and complexity of calculation, not only a strong mathematical and computer programming ability was required, but also a large number of long-term historical data needed to be used for calibrating the model parameters and improving its accuracy. In this study, three aspects were taken into consideration for the accuracy verification of the model. The first was to ensure the accuracy of the on-site monitoring data of PAHs and hydrology; the second was to calibrate the model parameters during the process of model construction; the third was to verify the accuracy of the simulation results after the calculation.

Collection and Experimental Analysis
In the summer of 2020, the investigation and sampling of hydrology, water quality and sediment in Pulandian Bay through a field survey were carried out. A total of 32 monitoring stations were set up. To verify the accuracy of the model simulation, three water-quality monitoring stations, three hydrological-monitoring stations (S3, S18 and S20), and two tidelevel observation stations (S11 and S15) were set up. The layout of the sampling stations is shown in Figure 1. A grab sampler (iron) was used to get the surface sediment (0-10 cm) samples, and a layered water sampler (stainless steel, 5 L) was used to get seawater samples. Hydrologic observations (flow velocity and direction) were collected every 10 min by using an Acoustic Doppler Current Profiler (Model: WHS600, depth unit size: 600 kHz).
Seawater samples were enriched by an envi-C18 solid phase membrane disk (flow rate: 12 mL/min−30 mL/min, vacuum filtration) extraction, wrapped in aluminum foil, sealed, and stored before testing (SPE samples). Sediment samples were freeze-dried (vacuum freeze dryer, 24 h) and ground (180 µm) after removal of impurities (e.g., plant residues and animal). Brown glass bottles were used for storage and then brought back for analysis. Sediment samples were weighed for about 20.0 g and then the following were added: anhydrous sodium sulfate 20.0 g, copper powder 20 g, 100 mL of extraction solution [acetone:n-hexane (V:V) = 1:1] and recovery indicators (D8-p, p'-DDT, CB155, CB65 and TBB). The samples were reflux-extracted through a Soxhlet apparatus for 24 h. The extract was evaporated to about 2 mL by rotary evaporation and stored. Water samples (SPE samples) were extracted with dichloromethane (3 times, 5 mL/time) and n-hexane (2 times, 5 mL/time), respectively, and the extract was concentrated to about 2.0 mL for further purification. Samples were purified by chromatography on silica gel/alumina columns. N-hexane (80 mL) and 7:3 v/v dichloromethane/n-hexane (100 mL) were used. The eluate was concentrated to 200 µL by high-purity nitrogen and spiked with terphenyl- (92.4%), respectively. All data had subtracted the blank. The limit of detection was set as 10 times the standard deviation of the spiked minimum-concentration blank (10 × SD). Assuming a concentration volume of 200 µL, the LOD was 0.02 ng/L for seawater and 0.05 ng/g (d.w.) for sediments.

Model Parameter Calibration
To improve the accuracy of the PAH numerical calculation model, the historical data (1999~2019) of seawater PAHs, which was monitored by the Dalian Marine Environment Monitoring Center Station of the State Oceanic Administration in Pulandian Bay, were used to calibrate the model parameters. The calibrated model parameters included the seawater viscosity coefficient (V e ), the transfer coefficient for the liquid or gas film (k l , k g ), the diffusion coefficient (D sw ), and the Chézy coefficient (Ψ).
Station S3 was used for model parameter calibration. The concentration of PAHs in seawater was taken as the judgment basis for the accuracy of model calibration results. According to the existing research [36], although the specific values of the above 5 parameters cannot be formulated directly, they are controlled by multiple factors such as topography and water depth. However, the value range of parameters can be estimated. Therefore, based on the existing research experience [37], the value range of each parameter was preliminarily set as V e ∈ [1.44, 1.57], D sw ∈ [0.12, 0.15], Ψ ∈ [0.02, 0.06], k g ∈ [0.73, 0.78], k l ∈ [4.12, 4.19].
In the process of parameter calibration, we found that the value of the seawater viscosity coefficient (V e ) had a strong impact on the concentration of PAHs in seawater simulated by the model [38]. Therefore, we encrypted the calibration of V e parameters (raised to 48 times). When the numerical simulation results of PAHs were basically consistent with their measured concentrations (in seawater), this indicated that the results of the model parameters were reasonable and the calibration work was completed. (See Figure 4 for further information). The parameter calibration results are shown in the Table 4 below.  To ensure the accuracy of the calibration results of model parameters, after determining the value of Ve, Dsw, Ψ, kg, and kl, the historical data (1999~2019) of seawater PAHs and hydrology data (including the tidal level, flow velocity and direction) were used to verify. The verification results showed that the predicted result of thePAHs numerical calculation model was in great agreement with the monitoring value, which proved that the model parameters were set appropriately ( Figure 5).

Model Accuracy Verification
To ensure the accuracy of the numerical simulation results, the numerical simulation model was verified based on the measured data (in 2020) of tide level, flow velocity, flow direction and ΣPAHs concentration in seawater continuously observed for 25 h at quality control stations (more information in Section 2.5.1 ). The results of model accuracy verification are shown in Figure 6. For further verification, the statistical method proposed by Willmott (1981) [39] was used to evaluate the model: where D i represented the observed value; − D represented the observed average value; i = 1, 2, . . ., n; Mi represented the simulated value; N represented the number of water samples; Skill represented the degree of correlation between the deviation of the observed value and observed average value and the deviation of the simulated value and the observed average value.
The model's performance was rated according to the Skill value: <0.2 was poor; 0.2-0.5 was medium; 0.5-0.7 was good; >0.7 was excellent. According to the verification results, the accuracy of the model was satisfactory: 0.83-0.88 for tide level; 0.81-0.84 for flow velocity; 0.85-0.89 for flow direction, and 0.77-0.83 for PAHs concentration in seawater indicated that the numerical hydrodynamic model was convincing and reliable.

Hydrodynamics Simulation Results
The hydrodynamics simulation results for Pulandian Bay are shown in Figure 7. The arrows denote the flow direction, and length is proportional to the flow velocity. The simulation results showed that the hydrodynamic conditions in Pulandian Bay were complex. From the bay-head to the bay-mouth, Pulandian Bay could be divided into two parts. With a high velocity, the bay-mouth area is relatively wide, and the flow direction is evenly and uniformly distributed; however, the bay-head area is narrow with a lower velocity, where the flow direction varies greatly and is distributed turbulently. In general, the hydrodynamics here had significant spatial distribution characteristics: the closer to the bay-head, the weaker the hydrodynamics.

The Simulation Results of PAHs
Time-series characteristic analysis: To explore the variation characteristics of ΣPAHs' concentration with time in Pulandian Bay, the model of PAHs was used to simulate the migration, transformation, and sedimentation process. The variation of PAHs simulated concentration with time at the bay-head (S20) and the bay-mouth (S3) was extracted from the model, and the concentration variation curve was drawn in Figure 8. From the curves, it could be learnt that with the simulated time going by, the concentration of ΣPAHs increased with fluctuations and finally reached a dynamic equilibrium (i.e., the concentration of PAHs fluctuates within a certain range). The concentration of ΣPAHs arrived at a stable level after 20 days (about 40 tidal cycles) at S3. However, the concentration of ΣPAHs did not reach dynamic equilibrium until 50 days (100 tidal cycles) at S20.
S3 was at the bay-mouth of Pulandian Bay, which had stronger hydrodynamics and ample water exchange. Consequently, the concentration of ΣPAHs at S3 was able to arrive at a dynamic equilibrium quickly, after only 20 days (40 tidal cycles). In contrast, the hydrodynamics of S20 was poor, located at the bay-head, where the ΣPAHs concentration kept increasing even after 40 days (80 tidal cycles). It seemed that this phenomenon was ordinary, but it was of great significance. The closer the offshore distance, the more complex were the hydrodynamic conditions and the weaker was the water-exchange capacity. Therefore, the diffusion of pollutants was more unfavorable. On the contrary, with a farther offshore distance, the hydrodynamic conditions would be stronger, which would be conducive to the rapid diffusion and dilution of pollutants.
It could be seen that the concentration of pollutants in the marine environment changed periodically, often due to tidal fluctuations. First-hand observation made it clear why it is difficult to get the real concentration of a pollutant in the marine environment through onetime monitoring, which was the largest problem for staff who had been engaged in marine environmental monitoring for a long time [40]. The migration and transformation process of pollutants in the marine environment is strongly affected by the hydrodynamics, which will be expanded further, especially in the bay or estuarine area. Time series monitoring is usually more convincing than spatial, but, due to the shortage of funds, manpower or equipment, it is very difficult and costly to implement regular monitoring at sea. Therefore, in order to explore the diffusion characteristic of pollutants and understand the tendencies of marine environmental quality, taking advantage of a numerical calculation model may be a more effective means.
Analysis of spatial distribution characteristics: After numerical simulation was started, the water in the calculation domain would exchange with the external to complete the PAHs transportation and migration. Specifically, with the tide rising, the external water (clean) flowed into the bay and mixed with the water init, resulting in the decreased concentration of ΣPAHs. After a certain time, the water exchange would be further completed towards the inner area of the bay. At ebb tide, the water with higher PAHs concentration in the bay would diffuse outside to complete the transportation and transfer. When the process of the next tidal cycle began, the cycle would be repeated. Figure 9 shows the concentration distribution of PAHs in high tide, falling tide, low tide, and rising tide during the numerical simulation, showing the spatial evolution pro-cess of PAHs migration and diffusion in Pulandian Bay. From the perspective of spatial distribution, the PAHs completed a water exchange from inside to outside under the force of hydrodynamics. The diffusion rate of PAHs varied significantly with the distance from the bay-head. Generally, the phenomenon showed a distribution trend of strengthening the hydrodynamic from the bay-head to the bay-Mouth. PAHs were able to stay in the bay-head area, where the hydrodynamic conditions were poor and the water exchange capacity was insufficient, that for an especially long time.

Analysis of PAHs in Seawater-Sediment Exchange
As described above, PAHs in the marine environment were easily adsorbed on suspended particulate materials and migrated with them. Therefore, when the suspended solids settled, the PAHs would also deposit and become part of the sediments. In order to explore the exchange characteristics of PAHs between seawater and sediment, based on the numerical simulation-calculating results, the distribution of seawater-sediment exchange rate of PAHs is illustrated ( Figure 10). Figure 10 shows that the sedimentation of PAHs was the primary process of seawatersediment exchange in Pulandian Bay, where the PAHs mainly migrated and transformed from seawater into sediment. Spatially, the results of the numerical simulations showed that the sedimentation rate exhibited an increasing distribution tendency from the bay-mouth (2 ng/m 2 ·day) to the bay-head (14 ng/m 2 ·day) in Pulandian Bay. According to the statistics of normal distribution, the sedimentation rate of PAHs (more than 60% grid units) was generally between 6-10 ng/m 2 ·day.
Compared with other bays (Guanabara Bay, Fraser River Basin e.g.,) in the world, the sedimentation rate of the PAHs was slightly high. After superposition, the total amount of PAHs transported from seawater to sediments every year was about 4.7 × 10 13 ng in Pulandian Bay. This might cause serious threats to aquaculture and human activities in nearshore waters. Because of the strong bioaccumulation function of PAHs, they might eventually cause damage to human health decades later. Therefore, in the following section, based on the monitoring data of PAHs, the numerical calculation model was coupled with the health risk assessment model so as to explore the extent to which pollution levels of PAHs in Pulandian Bay pose a risk to human health.

Health Risk Assessment of PAHs
As mentioned above, a single monitoring alongside survey results was not comprehensive enough, both temporally and spatially. However, most of the pollutant health risk assessment models were based on the identified and single survey/monitoring data results, which would inevitably introduce error and randomness. Therefore, to solve this difficulty, a dynamic model of health risk assessment for PAHs was built by coupling the health risk assessment model with the numerical calculation model. Depending on model parameter calibration and validation, it was possible to minimize the occasionality and uncertainty due to sampling and analysis processes, which could improve the authenticity and accuracy of evaluated conclusions.
As shown in Figure 11, the distribution of aquaculture in the Pulandian Bay was extremely intensive, especially in the nearshore waters, where the health risks of PAHs were so high that they could cause toxic effects on aquatic life or even pose a threat to human health through bioaccumulation. In this paper, it is suggested that the Marine Environmental Authority should pay more attention to control efforts to prevent PAH pollution in the head of Pulandian Bay. In order to minimize the pollution and impacts of PAHs on nearshore waters, monitoring amounts of PAHs pollutants or reference to the Pollutant Discharge Permit System might be an effective new measure [40].

Conclusions
The hydrodynamic conditions in the semi-closed, narrow Pulandian Bay were very complex and showed the conspicuous spatial distribution characteristics that the closer to the bay-head, the weaker the hydrodynamics.
The concentration of pollutants (PAHs) in the ocean fluctuated periodically with the tide, showing that hydrodynamic force had a significant impact on the diffusion of pollutants. The diffusion rate of pollutants in the area was faster with a stronger hydrodynamic force, and vice versa. This influence would be further amplified in the semi-closed narrow bay, which would lead to a lack of sufficient diffusion and dilution of pollutants in the bay-head area, resulting in a higher ecological risk.
The total amount of PAHs transported from seawater to sediments was 4.7 × 10 13 ng/year in Pulandian Bay. Arguably, it might cause serious threats to aquaculture and human activities in nearshore waters. Because of the strong bioaccumulation function of PAHs, the threats and damage may eventually be transferred to human health, decades later.